1 Introduction

This analysis looks at the field data collected for the following 2017 MLRA project: - MLRA 111B - Glynwood B-slope Erosion; Northeastern IN

Spatially disaggregate the existing SSURGO polygons for Glynwood B-slope map units using ArcSIE, in order to separate different soil erosion phases. The current SSURGO maps join issues at the SSA boundaries, due to different erosion phases. This project is deemed relevant due to the current interest in Soil Health. Distinguishing the difference in erosion phases may have minimal impact on the majority of soil interpretations, but is believed to be significant in distinguishing crop yields.

2 Data Preparation

library(aqp)
library(soilDB)

library(reshape2)
library(ggplot2)
library(gridExtra)
library(knitr)

library(cluster)
library(caret)
library(party)
library(vegan)

library(rgdal)
library(sp)
library(mapview)

2.1 Point Data

gp <- fetchNASIS()
s <- site(gp)[c("peiid", "pedon_id", "taxonname", "x", "y", "describer", "erocl")]
h <- horizons(gp)

vars <- c("peiid")
A_vars <- c("hzname", "hzdept", "clay", "texture", "m_hue", "m_value", "m_chroma")
Bt_vars <- c("hzdept", "clay", "texture", "m_value", "m_chroma")
carb_vars <- c("hzdept", "effervescence")
solum_vars <- c("hzdept")

h2 <- by(h, h[vars], function(x) data.frame(
  x[vars][1, ],
  x[A_vars][1, ],
  x[grepl("Bt", x$hzname), Bt_vars, ][1, ],
  x[x$effervescence %in% c("strong", "violent"), carb_vars][1,],
  x[grepl("^C|^2C|^3C", x$hzname), solum_vars][1]
  ))
h2 <- do.call("rbind", h2)
names(h2) <- c(vars, "hzname", "A_hzthk", "claytotest", "texture", "mxhue", "mxvalue", "mxchroma", "hzthk", "firstbtclay", "firstbttexture", "firstbtmxvalue", "firstbtmxchroma", "CaCO3Dp", "effervescence", "SolumDp")
names(s)[c(2:5, 7)]  <- c("upedonid", "soilname", "long", "lat", "EroClassFD")

gw <- merge(s, h2, by = "peiid", all.x = TRUE)
gw <- gw[order(gw$upedonid), ]
gw <- within(gw,{
             CaCO3Dp[is.na(CaCO3Dp)] <- 200
             SolumDp[is.na(SolumDp)] <- 200
             })

# projectmapunit data from NASIS
project <- get_projectmapunit_from_NASIS()
project_nodups <- project[!duplicated(project$nationalmusym), c("nationalmusym", "muname")]

# MUPOLYGONs for the Project
gw_pol <- readOGR(dsn = paste0(ownCloud, "glynwood.shp"), layer = "glynwood")

# Soil Survey Areas
ssa <- readOGR(dsn = paste0(geodata, "soils/soilsa_a_nrcs.shp"), layer = "soilsa_a_nrcs")

# Series Extent of Glynwood from SoilWeb
gw_series <- seriesExtent("Glynwood")

# field Data
gw2 <- read.csv(paste0(ownCloud, "Pts_gnbero_27Jan17.csv"))
vars <- c("upedonid", "EroClassSIE", "relpos", "SlopeSIE", "wetness", "PlanCrv", "ProfCrv", "maxcrv", "mincrv")
gw <- merge(gw, gw2[vars], by = "upedonid", all.x = TRUE)
gw <- gw[complete.cases(gw[c("lat", "long")]), ]

gw_sp <- gw
coordinates(gw_sp) <- ~ long + lat
proj4string(gw_sp) <- CRS("+init=epsg:4326")
gw_sp <- spTransform(gw_sp, CRS(proj4string(gw_pol)))

# spatial overlay field data with mupolygons and merge with nasis mapunits
vars <- c("AREASYMBOL", "nationalmu")
gw <- cbind(gw, over(gw_sp, gw_pol)[vars])
gw <- merge(gw, project_nodups, by.x = "nationalmu", by.y = "nationalmusym", all.x = TRUE, sort = FALSE)


# Extract erosion phases from NASIS and combine component and phase
ero_labels <- c("undisturbed", "slight", "moderate", "severe")

gw <- within(gw, {
  EroClassFD  = factor(EroClassFD, levels = 0:3, labels = ero_labels)
  EroClassSIE = factor(EroClassSIE, levels = 0:3, labels = ero_labels)
  EroClassFD2 = ifelse(EroClassFD == "severe", "severe", "slight")

  EroClassNASIS = NA
  EroClassNASIS[grepl("eroded", muname)] = "eroded"
  EroClassNASIS[grepl("sev.|severely", muname)] = "sev.eroded"
  EroClassNASIS[!grepl("eroded", muname)] = "non.eroded"
  
  soilname2 = soilname
  soilname2 = ifelse(soilname2 %in% c("Glynwood", "Morley", "Shinrock", "Rawson", "Mississinewa"), "Glynwood", soilname2)
  soilname2 = ifelse(soilname2 %in% c("Blount", "Elliott"), "Blount", soilname2)
  soilname2 = ifelse(soilname2 %in% c("Pewamo", "Pandora", "Mermill"), "Pewamo", soilname2)
  soilname3 = paste0(soilname2, ifelse(soilname2 == "Glynwood", paste0("-", EroClassFD2), ""))
  })
gw <- transform(gw,
                rgb = munsell2rgb(mxhue, mxvalue, mxchroma, return_triplets = TRUE)
                )
gw_sp <- gw
coordinates(gw_sp) <- ~ long + lat
proj4string(gw_sp) <- CRS("+init=epsg:4326")
gw_sp <- spTransform(gw_sp, CRS(proj4string(gw_pol)))

2.2 Spatial Data

The geodata from the Glynwood points was extracted from several rasters at various resolutions. The data using to generate the ArcSIE model came from a DEM with a resolution of 15-feet. The other used came from the 10-meter USGS NED, which was primarily resampled from LiDAR.

# Extract data from rasters
library(raster)

# NW files
fd <- paste0(geodata, "project_data/11FIN/PointDataEval/")
dd <- c("slope10",
        "procrv10",
        "plncrv10",
        "maxcrv10",
        "mincrv10",
        "relpos_r5",
        "wetness_mp"
        )
fp <- paste0(fd, "Mosaic_NW_pts/Derivatives/", dd, "/", "w001001.adf")
rs <- stack(fp)
names(rs) <- dd
proj4string(rs) <- CRS("+init=epsg:2965")
gd_nw <- extract(rs, gw_sp, df = TRUE, sp = TRUE)@data
gd_nw <- subset(gd_nw, !is.na(slope10))

# SE files
fp <- paste0(fd, "Mosaic_SE_pts/Derivatives/", dd, "/", "w001001.adf")
rs <- stack(fp)
names(rs) <- dd
proj4string(rs) <- CRS("+init=epsg:2965")
gd_se <- extract(rs, gw_sp, df = TRUE, sp = TRUE)@data
gd_se <- subset(gd_se, !is.na(slope10))

gd15ft <- rbind(gd_nw, gd_se)
rm(gd_nw, gd_se)
write.csv(gd15ft, file = paste0(ownCloud, "geodata_15ft_extract.csv"))


# Region 11 files

subset_rasters <- function(input, output) {
  cat(paste0(input, "\n"))
  gdal_translate(
    src_dataset = input,
    dst_dataset = output,
    projwin = c(bb[1,1], bb[2,2], bb[1,2], bb[2,1]),
    of = "GTiff",
    a_nodata = -99999,
    overwrite = TRUE,
    verbose = TRUE
    )
  }

warp_rasters <- function(input, output){
  cat(paste0(input,"\n"))
  gdalwarp(
    srcfile = input,
    dstfile = output,
    te = bb,
    s_srs = CRSargs(CRS("+init=epsg:5070")),
    t_srs = CRSargs(CRS("+init=epsg:5070")),
    r = "bilinear",
    tr = c(10, 10),
    of = "GTiff",
    overwrite = TRUE,
    verbose = TRUE 
    )
  }

fd <- paste0(geodata, "project_data/11FIN/sdat/")
dd2 <- c("ned10m_11FIN.sdat",
        "ned10m_11FIN_aspect5.sdat",
        "ned10m_11FIN_slope5.sdat",
        "ned10m_11FIN_cupro5.sdat",
        "ned10m_11FIN_cutan5.sdat",
        "ned30m_11FIN_mvalleys.sdat",
        "ned30m_11FIN_wetness.sdat",
        "ned30m_11FIN_z2stream.sdat",
        "nlcd30m_11FIN_lulc2011.sdat"
        )
dd_names <- c("elev", "aspect5", "slope5", "kp5", "kt5", "mvalley", "wetness2", "z2streams", "lulc")
dd <- paste0(fd, dd2)

input <- dd
output <- paste0(geodata, "project_data/11FIN/glynwood/", gsub(".sdat", ".tif", dd2))
bb <- bbox(gw_sp <- spTransform(gw_sp, CRS("+init=epsg:5070")))

mapply(subset_rasters, input, output)
mapply(subset_rasters, input = paste0(geodata, "project_data/11FIN/nlcd30m_11FIN_lulc2011.tif"), output = paste0(geodata, "project_data/11FIN/glynwood/nlcd30m_11FIN_lulc2011.tif"))

mapply(warp_rasters, input = output, output = gsub(".tif", "2.tif", gsub("30m", "10m", output)))

dd <- output
 
rs10m <- stack(dd[grepl("10m", dd)])
names(rs10m) <- dd_names[1:5]
rs30m <- stack(dd[grepl("30m", dd)])
names(rs30m) <- dd_names[6:9]

rs10m <- stack(gsub(".tif", "2.tif", gsub("30m", "10m", output)))
names(rs10m) <- dd_names

gd10m <- as.data.frame(extract(rs10m, gw_sp, df = TRUE, sp = TRUE))
gd30m <- as.data.frame(extract(rs30m, gw_sp, df = TRUE))
gw <- cbind(gd10m, gd30m[, -1])
rm(gd10m, gd30m)


# Save data
save(gw, gw_sp, gw_pol, gw_series, ssa, ero_labels, file = paste0(ownCloud, "Pts_gnbero_27Jan17_geodata2.RData"))

3 Map Units vs Field Observations

# Load cached dataset
load(paste0(ownCloud, "Pts_gnbero_27Jan17_geodata2.RData"))

ssa <- subset(ssa, areasymbol %in% unique(gw$AREASYMBOL))
ssa <- spTransform(ssa, CRS("+init=epsg:4326"))

# Create interactive map
mapView(gw_series) + ssa + gw_sp

vals2 <- c("EroClassFD", "EroClassNASIS", "nationalmu", "AREASYMBOL")
gw_sub <- gw[vals2]

# Frequency of field observation vs map unit
# Duplicate the data for each REASYMBOL and relabel MLRA
gw_sub2 <- by(gw_sub, gw_sub$nationalmu, function(x) {
  x[vals2][1, ]
  x[, "AREASYMBOL"] <- "MLRA"
  return(x)
  })
gw_sub2 <- do.call("rbind", gw_sub2)
gw_sub <- rbind(gw_sub, gw_sub2)
gw_sub$natmuSsaEro <-  with(gw_sub,
                            paste0(nationalmu, "-", AREASYMBOL, "-", EroClassNASIS)
                            )
test <- xtabs(~ natmuSsaEro + EroClassFD, data = gw_sub)
kable(test, caption = "Frequence by MUSYM-SSA-EROSION")

Frequence by MUSYM-SSA-EROSION
	undisturbed	slight	moderate	severe
2t6ll-IN009-sev.eroded	13	14	17	3
2t6ll-IN053-sev.eroded	4	11	17	11
2t6ll-IN075-sev.eroded	6	8	4	5
2t6ll-IN179-sev.eroded	1	9	6	10
2t6ll-MLRA-sev.eroded	24	42	44	29
2t6lm-IN009-sev.eroded	2	11	10	1
2t6lm-IN053-sev.eroded	2	4	8	11
2t6lm-IN075-sev.eroded	3	1	11	9
2t6lm-IN179-sev.eroded	9	8	0	13
2t6lm-MLRA-sev.eroded	16	24	29	34
2v4bn-IN069-eroded	4	3	6	4
2v4bn-IN179-eroded	0	3	4	6
2v4bn-MLRA-eroded	4	6	10	10
2v4bp-IN179-eroded	0	3	0	2
2v4bp-MLRA-eroded	0	3	0	2
5jjt-IN035-sev.eroded	1	1	2	9
5jjt-MLRA-sev.eroded	1	1	2	9
NA-NA-non.eroded	6	9	10	34

kable(round(prop.table(test, 1) * 100), caption = "Percent by MUSYM-SSA-EROSION")

Percent by MUSYM-SSA-EROSION
	undisturbed	slight	moderate	severe
2t6ll-IN009-sev.eroded	28	30	36	6
2t6ll-IN053-sev.eroded	9	26	40	26
2t6ll-IN075-sev.eroded	26	35	17	22
2t6ll-IN179-sev.eroded	4	35	23	38
2t6ll-MLRA-sev.eroded	17	30	32	21
2t6lm-IN009-sev.eroded	8	46	42	4
2t6lm-IN053-sev.eroded	8	16	32	44
2t6lm-IN075-sev.eroded	12	4	46	38
2t6lm-IN179-sev.eroded	30	27	0	43
2t6lm-MLRA-sev.eroded	16	23	28	33
2v4bn-IN069-eroded	24	18	35	24
2v4bn-IN179-eroded	0	23	31	46
2v4bn-MLRA-eroded	13	20	33	33
2v4bp-IN179-eroded	0	60	0	40
2v4bp-MLRA-eroded	0	60	0	40
5jjt-IN035-sev.eroded	8	8	15	69
5jjt-MLRA-sev.eroded	8	8	15	69
NA-NA-non.eroded	10	15	17	58

Several of counties phased severely eroded, are not dominanted by field observations classified as severely eroded.

4 Accuracy Assessment of the ArcSIE Predictions

cm <- confusionMatrix(data = gw$EroClassSIE, reference = gw$EroClassFD)
print(cm)

## Confusion Matrix and Statistics
## 
##              Reference
## Prediction    undisturbed slight moderate severe
##   undisturbed           0      0        0      0
##   slight                9      6       10      9
##   moderate             27     34       51     37
##   severe               12     31       21     45
## 
## Overall Statistics
##                                          
##                Accuracy : 0.3493         
##                  95% CI : (0.2947, 0.407)
##     No Information Rate : 0.3116         
##     P-Value [Acc > NIR] : 0.09322        
##                                          
##                   Kappa : 0.0862         
##  Mcnemar's Test P-Value : 1.115e-14      
## 
## Statistics by Class:
## 
##                      Class: undisturbed Class: slight Class: moderate
## Sensitivity                      0.0000       0.08451          0.6220
## Specificity                      1.0000       0.87330          0.5333
## Pos Pred Value                      NaN       0.17647          0.3423
## Neg Pred Value                   0.8356       0.74806          0.7832
## Prevalence                       0.1644       0.24315          0.2808
## Detection Rate                   0.0000       0.02055          0.1747
## Detection Prevalence             0.0000       0.11644          0.5103
## Balanced Accuracy                0.5000       0.47891          0.5776
##                      Class: severe
## Sensitivity                 0.4945
## Specificity                 0.6816
## Pos Pred Value              0.4128
## Neg Pred Value              0.7486
## Prevalence                  0.3116
## Detection Rate              0.1541
## Detection Prevalence        0.3733
## Balanced Accuracy           0.5880

test <- as.data.frame(cm$table)

ggplot(test, aes(x = Reference, y = Freq, fill = Prediction)) +
  geom_bar(stat = "identity") +
  coord_flip()

The accuracy of the current ArcSIE model appears to be low, according to several metrics. The positive predictive value for the severe class is < 50%.

4.1 Boxplots

soil_vals <- c("hzthk", "SolumDp", "CaCO3Dp", "claytotest", "firstbtclay", "mxvalue", "mxchroma")
geo_vals1 <- c("SlopeSIE", "ProfCrv", "PlanCrv", "relpos", "wetness")
geo_vals2 <- c("slope5", "kt5", "kp5", "z2streams", "wetness2", "mvalley")

vals <- c(soil_vals, geo_vals1, geo_vals2)
gw <- gw[complete.cases(gw[c("EroClassFD", soil_vals)]), ]

gw_lo1 <- melt(gw, id.vars = "EroClassFD", measure.vars = vals)
gw_lo2 <- melt(gw, id.vars = "EroClassSIE", measure.vars = vals)

names(gw_lo1)[1] <- "EroClass"
gw_lo1$method <- "FD"
names(gw_lo2)[1] <- "EroClass"
gw_lo2$method <- "SIE"
gw_lo <- rbind(gw_lo1, gw_lo2)
gw_lo <- subset(gw_lo, !is.na(EroClass))

gw_lo <- na.exclude(gw_lo)
ggplot(gw_lo, aes(x = EroClass, y = value)) +
  geom_boxplot() +
  facet_wrap(~ paste(variable, method), scales="free", ncol = 4) +
  coord_flip()

An exploratory analysis shows a considerable amount of overlap exists between the field determined (FD) erosion classes and measurable soil properties. In comparison the FD and SIE (Soil Inference Engine) erosion classes show different patterns within the boxplots, further suggesting that the SIE classes aren’t capturing the field observations accurately. The most important feature to highlight is that the trends between the SIE classes and digital elevation model (DEM) derivatives (i.e. slope) don’t match those observed for the FD classes. This mismatch suggests that the membership functions for the SIE classes are a poor fit, and should be redefined to more accurately represent the relationship between the FD classes and DEM derivatives.

4.2 Scatterplots

soil_vals2 <- c("hzthk", "SolumDp", "CaCO3Dp", "claytotest", "firstbtclay") # excluded color, only observed a narrow range thus small differences swamp everthing else
vals <- c(soil_vals2)

test <- gw[, vals]
test_d <- daisy(scale(test), metric = "gower")
test_mds <- metaMDS(test_d, distance = "gower", autotransform = FALSE, trace = FALSE)
test_pts <- cbind(as.data.frame(test_mds$points), EroClassFD = gw$EroClassFD)

g1 <- ggplot(gw, aes(x = hzthk, y = SolumDp, color = EroClassFD)) +
  geom_point(cex = 2, alpha = 0.75) +
  theme(aspect.ratio = 1)
g2 <- ggplot(test_pts, aes(x = MDS1, y = MDS2, color = EroClassFD)) +
  geom_point(cex = 2, alpha = 0.75) +
  theme(aspect.ratio = 1)
grid.arrange(g1, g2, ncol = 2)

According to the scatterplot above it appears that only the severe and slight classes are separatable. The moderate erosion class seems to overlap the most with slight. The overlap in the FD classes is likely due to bias within and between the soil scientists who collected the data. Both the 15-feet and 10-meter DEM derivatives were evaluated, but the results are similar.

4.3 Classification Tree

test <- subset(gw, !is.na(EroClassFD))

test_ct <- ctree(EroClassFD ~ ., data = test[, c("EroClassFD", soil_vals)])
plot(test_ct)

cm <- confusionMatrix(data = predict(test_ct, type = "response"), reference = test$EroClassFD)
print(cm)

## Confusion Matrix and Statistics
## 
##              Reference
## Prediction    undisturbed slight moderate severe
##   undisturbed          37      5        1      0
##   slight                8     62        3      1
##   moderate              5     14       66      4
##   severe                0      2       23    104
## 
## Overall Statistics
##                                           
##                Accuracy : 0.803           
##                  95% CI : (0.7563, 0.8442)
##     No Information Rate : 0.3254          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.729           
##  Mcnemar's Test P-Value : NA              
## 
## Statistics by Class:
## 
##                      Class: undisturbed Class: slight Class: moderate
## Sensitivity                      0.7400        0.7470          0.7097
## Specificity                      0.9789        0.9524          0.9050
## Pos Pred Value                   0.8605        0.8378          0.7416
## Neg Pred Value                   0.9555        0.9195          0.8902
## Prevalence                       0.1493        0.2478          0.2776
## Detection Rate                   0.1104        0.1851          0.1970
## Detection Prevalence             0.1284        0.2209          0.2657
## Balanced Accuracy                0.8595        0.8497          0.8073
##                      Class: severe
## Sensitivity                 0.9541
## Specificity                 0.8894
## Pos Pred Value              0.8062
## Neg Pred Value              0.9757
## Prevalence                  0.3254
## Detection Rate              0.3104
## Detection Prevalence        0.3851
## Balanced Accuracy           0.9218

An analysis of the EroClassFD above with a classification tree is an attempt to discern the hierachical structuce within the data. The results show Ap thickness (hzthk) and clay content (claytotest) are the first splits. The trees structure follows the logic described in the erosion indicators guide developed for this project. The overall accuracy for the tree is 0.8.

5 Hierachical Clusters

In order to see if more separation can be achieved amongst the erosion classes a hierachical classifition was peformed. Four unsupervised classes were chosen and manually matched to the FD classes.

test_c <- hclust(test_d, method = "ward")
plot(test_c, labels = gw$upedonid)
rect.hclust(test_c, k = 4)

clusters <- cbind(gw, 
                  test_pts[, 1:2], 
                  clusters = factor(cutree(test_c, k = 4),
                                    levels = c(2, 3, 1, 4),
                                    labels = ero_labels
                                    )
                  )
clusters <- cbind(gw, 
                  test_pts[, 1:2], 
                  clusters = factor(cutree(test_c, k = 4),
                                    levels = c(4, 1, 2, 3),
                                    labels = ero_labels
                                    )
                  )

xtabs(~ EroClassFD + clusters, data = clusters)

##              clusters
## EroClassFD    undisturbed slight moderate severe
##   undisturbed          31     14        5      0
##   slight               11     54       18      0
##   moderate              5     28       43     17
##   severe                0      5       34     70

5.1 Scatter Plots

g1 <- ggplot(clusters, aes(x = MDS1, y = MDS2, col = EroClassFD)) +
  geom_point(cex = 2, alpha = 0.75) +
  theme(aspect.ratio = 1)
g2 <- ggplot(clusters, aes(x = MDS1, y = MDS2, col = clusters), main = "test") +
  geom_point(cex = 2, alpha = 0.75) +
  theme(aspect.ratio = 1)
grid.arrange(g1, g2, ncol = 2)

In comparison the hierarchical clusters have less overlap when viewed along the multidimensional scaled axes, but still does not seem to separate the moderate class.

5.2 Box Plots

gw_lo3 <- melt(clusters, id.vars = "clusters", measure.vars = c(soil_vals, geo_vals2))

names(gw_lo3)[1] <- "EroClass"
gw_lo3$method <- "clusters"
gw_lo1 <- subset(gw_lo1, ! variable %in%  c("relpos", "wetness", "SlopeSIE"))
gw_lo <- rbind(gw_lo1, gw_lo3)

ggplot(gw_lo, aes(x = EroClass, y = value)) +
  geom_boxplot() +
  facet_wrap(~ paste(variable, method), scales="free", ncol = 4) + 
  coord_flip()

A comparison of the FD and cluster classes shows that the clusters do a good job replicating the patterns found in the boxplots.

5.3 Classification Tree

test2 <- ctree(clusters ~ ., data = clusters[, c("clusters", soil_vals)])
plot(test2)

confusionMatrix(data = predict(test2, type = "response"), reference = clusters$clusters)

## Confusion Matrix and Statistics
## 
##              Reference
## Prediction    undisturbed slight moderate severe
##   undisturbed          47      0        0      0
##   slight                0     91        3      1
##   moderate              0      7       89      1
##   severe                0      3        8     85
## 
## Overall Statistics
##                                          
##                Accuracy : 0.9313         
##                  95% CI : (0.8988, 0.956)
##     No Information Rate : 0.3015         
##     P-Value [Acc > NIR] : < 2.2e-16      
##                                          
##                   Kappa : 0.9065         
##  Mcnemar's Test P-Value : NA             
## 
## Statistics by Class:
## 
##                      Class: undisturbed Class: slight Class: moderate
## Sensitivity                      1.0000        0.9010          0.8900
## Specificity                      1.0000        0.9829          0.9660
## Pos Pred Value                   1.0000        0.9579          0.9175
## Neg Pred Value                   1.0000        0.9583          0.9538
## Prevalence                       0.1403        0.3015          0.2985
## Detection Rate                   0.1403        0.2716          0.2657
## Detection Prevalence             0.1403        0.2836          0.2896
## Balanced Accuracy                1.0000        0.9419          0.9280
##                      Class: severe
## Sensitivity                 0.9770
## Specificity                 0.9556
## Pos Pred Value              0.8854
## Neg Pred Value              0.9916
## Prevalence                  0.2597
## Detection Rate              0.2537
## Detection Prevalence        0.2866
## Balanced Accuracy           0.9663

In comparision, the classification tree for the clusters splits primarily on the CaCO3 and solum depths, presumable due to the narrow range in Ap thickness.

6 Statistical Modeling

Below several statistical models were evaluated to see if a more accurate model could be developed.

6.1 FD Classes vs DEM Derivatives

test3 <- ctree(EroClassFD ~ ., data = gw[, c("EroClassFD", geo_vals2)])
plot(test3)

cm_ct <- confusionMatrix(data = predict(test3, type = "response"), reference = gw$EroClassFD)
round(cm_ct$overall, 2)

##       Accuracy          Kappa  AccuracyLower  AccuracyUpper   AccuracyNull 
##           0.38           0.11           0.32           0.43           0.33 
## AccuracyPValue  McnemarPValue 
##           0.03           0.00

test3 <- cforest(as.factor(EroClassFD) ~ ., data = gw[, c("EroClassFD", geo_vals2)])
varimp(test3)

##      slope5         kt5         kp5   z2streams    wetness2     mvalley 
## 0.037024390 0.019056911 0.005479675 0.001138211 0.016569106 0.000699187

cm_cf <-confusionMatrix(data = predict(test3, type = "response", OOB = TRUE), reference = gw$EroClassFD)
round(cm_cf$overall, 2)

##       Accuracy          Kappa  AccuracyLower  AccuracyUpper   AccuracyNull 
##           0.37           0.14           0.32           0.43           0.33 
## AccuracyPValue  McnemarPValue 
##           0.04           0.03

Neither a classification tree or forest were capiable of achieving a higher accuracy than the SIE model.

6.2 Clusters vs DEM Derivatives

test4 <- ctree(clusters ~ ., data = clusters[, c("clusters", geo_vals2)])
plot(test4)

cm_ct <- confusionMatrix(data = predict(test4, type = "response"), reference = clusters$clusters)
round(cm_ct$overall, 2)

##       Accuracy          Kappa  AccuracyLower  AccuracyUpper   AccuracyNull 
##           0.38           0.12           0.32           0.43           0.30 
## AccuracyPValue  McnemarPValue 
##           0.00           0.00

test4 <- cforest(clusters ~ ., data = clusters[, c("clusters", geo_vals2)])
varimp(test4)

##       slope5          kt5          kp5    z2streams     wetness2 
## 0.0186016260 0.0343252033 0.0004227642 0.0017073171 0.0084390244 
##      mvalley 
## 0.0131544715

cm_cf <- confusionMatrix(data = predict(test4, type = "response", OOB=TRUE), reference = clusters$clusters)
round(cm_cf$overall, 2)

##       Accuracy          Kappa  AccuracyLower  AccuracyUpper   AccuracyNull 
##           0.36           0.12           0.31           0.42           0.30 
## AccuracyPValue  McnemarPValue 
##           0.01           0.39